# -----------Masonry infill -------------------------------------------------------------------------
# -- Gianrocco Mucedero PhD Student------------------------------------------------------------------
# --Understanding and Managing Extremes--------------------------------------------------------------
# -----------IUSS Pavia -----------------------------------------------------------------------------
# -----------23-09-2019 -----------------------------------------------------------------------------
# ---------------------------------------------------------------------------------------------------
# -- Description of the Parameters -------
# eleTag 	label of panel
# typ		type of infill model (single/double/triple)
# nds 		list of nodes [list TL TR BR BL] (i.e always start numbering from top left corner)
# B 		Bay width (centrelines) [mm]
# H 		Height (centrelines) [mm]
# hc 		Column height [mm]
# bc 		Column width [mm]
# bb        Beam width [mm]
# hb        Beam height [mm]
# hp        Opening heigh  [mm]
# lp        Opening length [mm]
# Ec 		Concrete Elastic Modulus [MPa]
# tw 		Wall thickness [mm]
# Ewv 		Vertical Secant Modulus [MPa]
# Ewh 		Horizontal secant modulus [MPa]
# Gw 		Shear Modulus [MPa]
# v 		Poissons Ratio
# fwv 		Vertical compressive strength [MPa]
# fwu 		Sliding shear resistance of mortar joints [MPa]
# fws 		Shear resistance under diagonal compression [MPa]
# flat	    Lateral compression strength of masonry	[MPa]
# sig_v		Vertical Compression due to gravity loading [MPa]
# sigm_cr   Cracking stregth of masonty [MPa]
# --------------------------------------------------------------------------------
# New inputs in addition to the original ones (Gerard O'Reilly function):
# --------------------------------------------------------------------------------
# bb        Beam width [mm]
# hb        Beam height [mm]
# hp        Opening heigh  [mm]
# lp        Opening length [mm]
# flat	    Lateral compression strength of masonry	[MPa]
#put flag option to Set: 
#        WIDTH OF STRUT: $StrutWidthflag==
#                        1 Bertoldi 
#                        2 Paulay and Priestley 1992
#                        3 Holmes 1961 
#                        4 Liauw and Kwan 1984
#                        5 Mainstone 1974 (Adopted in FEMA 306)
#                        6 Stafford  Smith 1961 : there are 2 different options (one is equal to case 2 Paulay and Priestley 1992) and the other is 0.1*dw 
#                        7 Decanini-Fantin uncracked masonry infill
#                        8 Decanini-Fantin cracked masonry infill
#                        9 Cavalieri et al. 2015
#put flag option to Set: 
#        REDUCTION COEFFICIENT STIFFNESS: $Reductionflag==
#                        0 No reduction 
#                        1 Dawe and Seah 
#                        2 Imai and Miyamoto 
#                        3 Tasnini and Mohebkhan 
#                        4 Decanini et al.
#                        5 Asteris et al.
#put flag option to Set: 
#        FAILURE MECHANISM: $Criticalstressflag==
#                        1 Decanini Fantini-Bertoldi et al. 1987
#                        2 Decanini Fantini-Bertoldi et al. 1987 with a reduction of fwu (shear sliding of bed joint ) according to Morandi et al. 2018 and 2018 and EC6
#                        3 Paulay Priestley and Calvi 
#                        4 Paulay Priestley and Calvi  with a reduction of fwu (shear sliding of bed joint ) according to Morandi et al. 2018 and 2018 and EC6
#						 5 FEMA 306
#                        6 Eurocode EC8-Part 1/EC6 but without Fvo REDUCTION
#                        7 Eurocode EC8-Part 1/EC6 with Fvo REDUCTION
#put flag option to Set: 
#        BACKBONE HYSTERESIS: $Backboneflag==
#                        1 Bertoldi
#                        2 Panagiotakos and Fardis
#                        3 De Risi et al. 
#                        4 Bertoldi + Sassun et. al 
#                       
#put flag option to Set: 
#        HYSTERETIC BEHAVIOUR OF MASONRY INFILL: $Hysteresticflag==
#                        1 Coefficients for Unloading/Reloading Stiffness degradation and Coefficients for Strength degradation ==0
#                        2 Coefficients for Unloading/Reloading Stiffness degradation and Coefficients for Strength degradation based on N. Mohammad Noh et al. / Engineering Structures 150 (2017) 599–621

# --------------------------------------------------------------------------------
# Start procedure to define Backbone hysteretic behaviour of masonry infill 
# --------------------------------------------------------------------------------
# Create a function that creates a Pinching4 material in a single line

proc infill {eleTag typ nds B H hc bc hb bb hp lp Ec MasonryInfillType StrutWidthflag Reductionflag Criticalstressflag Backboneflag Hysteresticflag pfile} {
# Flag for : Reduction due to opening, Strut width relation, Failure mechanism, Backbone model, Hysteretic behaviour.
source Units.tcl
if {$MasonryInfillType==1} {
# Weak Infill (Parameters from Hak et al. 2012)
set tw		80.00
set Ewv	    1873.00
set Ewh	    991.00
set Gw	    1089.00
set v       0.250
set fwv	    2.020
set fwu	    0.44
set fws 	0.55
set flat	1.18
set sig_v   0.0
set sigm_cr 0.0
}
if {$MasonryInfillType==2} {
# Medium Infill (Parameters from Hak et al. 2012)
set tw		240.00
set Ewv 	1873.00
set Ewh 	991.00
set Gw		1089.00
set v       0.250
set fwv  	1.500
set fwu 	0.250
set fws 	0.310
set flat	1.110
set sig_v   0.0
set sigm_cr 0.0
}
if {$MasonryInfillType==3} {
# Strong Infill (Parameters from Hak et al. 2012)
set tw		300.00
set Ewv 	3240.00
set Ewh 	1050.00
set Gw  	1296.00
set v       0.250
set fwv 	3.510
set flat	1.500
set fwu 	0.300
set fws 	0.360
set sig_v   0.0
set sigm_cr 0.0
}
if {$MasonryInfillType==4} {
# Morandi et al. 2018 
set tw		350.00
set Ewv 	5299.00
set Ewh 	494.00
set Gw  	2120.00
set v       0.25
set fwv 	4.640
set flat	1.080
set fwu 	0.359
set fws 	0.0
set sig_v   0.0
set sigm_cr 0.0
}
if {$MasonryInfillType==5} {
# Cavalieri et al. 2014 Clay masonry S1B-S2B
set tw	    150.00	
set Ewv 	6401.00
set Ewh 	5038.00
set Gw  	2547.00
set v       0.25
set fwv 	8.660
set flat	4.180
set fwu 	1.07
set fws 	0.0
set sig_v   0.0
set sigm_cr 0.0
}

if {$MasonryInfillType==6} {
# Cavalieri et al. 2014 Calcarenite masonry S2-A
set tw	    210.00	
set Ewv 	7106.00
set Ewh 	9528.00
set Gw  	2937.00
set v       0.250
set fwv 	4.570
set flat	3.920
set fwu 	0.89
set fws 	0.0
set sig_v   0.0
set sigm_cr 0.0
}

set pi 3.14
# --------------------------------------
# Additional TERMS detetrmination
# --------------------------------------
set Ac      [expr $hc*$bc];           # mm2 , area of column 
set Ab      [expr $hb*$bb];           # mm2 , area of beam
set Ib 		[expr $bb*pow($hb,3)/12]; # mm4 , inertia of the beam 
set Ic 		[expr $bc*pow($hc,3)/12]; # mm4 , inertia of the column
set Bw 		[expr $B-$hc];            # mm
set Hw		[expr $H-$hb];            # mm
set dw		[expr sqrt($Bw*$Bw+$Hw*$Hw)]; # mm, length of the strut 
set theta 	[expr atan2($Hw,$Bw)];    # rad
set theta2   [expr $theta*180/$pi];
set Ewtheta	[expr 1/(pow(cos($theta),4)/$Ewh+pow(sin($theta),4)/$Ewv+pow(cos($theta),2)*pow(sin($theta),2)*(1.0/$Gw-2.0*$v/$Ewh))];
set lambdaH	[expr  $H*(pow(($Ewtheta*$tw*sin(2*$theta)/4/$Ec/$Ic/$Hw),0.25))];
set z		[expr  $pi/2/$lambdaH*$H]; # mm
set s2		[expr  $z/3]; # in mm
set s3		[expr  $z/2]; # in mm
puts [format "Area element mm2  : 
Ac=%.3f Ab=%.3f  " $Ac $Ab]
puts [format "Inerzia element mm4  : 
Ic=%.3f Ib=%.3f  " $Ic $Ib]
puts [format "Infill properties mm :
Height= %.3f  Length= %.3f" $Hw $Bw ]
puts [format "Diagonal length mm :  
dw= %.3f" $dw $theta2 ]
puts [format "Angle of strut: 
Theta = %.3f " $theta2 ]
puts [format "Elastic modulus in diagonal direction MPa: 
Ewtheta =%.3f" $Ewtheta ]
puts [format "Lamda H: 
lambdaH= %.3f" $lambdaH];
puts [format "Contact length mm : 
Z= %.3f" $z];
puts [format "Strut off diagonal position 2 strut model mm : 
s2= %.3f" $s2];
puts [format "Strut off diagonal position 3 strut model mm : 
s3= %.3f" $s3];

# ################# Next step to be implemented ##############################
set lamdbacolumn [expr  $H*(pow(($Ewtheta*$tw*sin(2*$theta)/4/$Ec/$Ic/$Hw),0.25))];  #Condisering column for relative stiffnees panel-frame
set lamdbabeam   [expr  $H*(pow(($Ewtheta*$tw*sin(2*$theta)/4/$Ec/$Ib/$Hw),0.25))];  #Condisering beam for relative stiffnees panel-frame
set zcolumn      [expr  $pi/2/$lamdbacolumn]; # mm   Contact length column-panel 
set zbeam        [expr  $pi/2/$lamdbabeam];   # mm   Contact length beam-panel 
# ############################################################################


# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# DETERMINE THE REDUCTION STIFFNESS COEFFICIENT 
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
#put flag option to Set: 
#        Reduction stiffness coefficient: 
#                        0 No reduction 
#                        1 Dawe and Seah 
#                        2 Imai and Miyamoto 
#                        3 Tasnini and Mohebkhan 
#                        4 Decanini et al.
#                        5 Asteris et al.

if {$Reductionflag >0} { 
set alphaA [expr  100.00*($lp*$hp)/($Bw*$Hw)];
set alphaL [expr  100.00*$lp/$Bw];
set alphaH [expr  100.00*$hp/$Hw];

puts [format "Reduction opening: alphaA=%.3f alphaL= %.3f  alphaH= %.3f"    $alphaA $alphaL $alphaH];
}
if {$Reductionflag==0} {
set rp [expr (1.0)];

puts [format "Reduction opening-No opening: 
Rp=%.3f" $rp]
}
# 1 Dawe and Seah
if {$Reductionflag==1} {
set rp [expr  1-(1.500 * ($alphaL/100))];

puts [format "Reduction opening  Dawe and Seah: 
Rp=%.3f" $rp]
}
# 2 Imai and Miyamoto
if {$Reductionflag==2} {
set par1 [expr 1-(0.01 * $alphaL)];
set par2 [expr 1-(0.1*pow($alphaA,0.5))];
set rp  [expr min($par1,$par2)];

puts [format "Reduction opening Imai and Miyamoto : 
Rp=%.3f" $rp]
}

# 3 Tasnini and Mohebkhan
if {$Reductionflag==3} {
set rp [expr 1-2.238*($alphaA/100)+1.49*pow(($alphaA/100),2)];

puts [format "Reduction opening Tasnini and Mohebkhan : 
Rp=%.3f" $rp]
}
# 4 Decanini et al.
if {$Reductionflag==4} { 
set e  2.7182818284590
set rp  [expr 0.55*pow($e,(-0.035*$alphaA))+0.44*pow($e,(-0.025*$alphaL))];

puts [format "Reduction opening Decanini et al. : 
Rp=%.3f" $rp]
}
# 5 Asteris et al.
if {$Reductionflag==5} {
set rp  [expr 1-2*(pow(($alphaA/100),0.54))+ pow(($alphaA/100),1.14)];

puts [format "Reduction opening Asteris et al. : 
Rp=%.3f" $rp]
}
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# DETERMINE WIDTH OF STRUT
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
#put flag option to Set: 
#        WIDTH OF STRUT: 
#                        1 Bertoldi 
#                        2 Paulay and Priestley 1992
#                        3 Holmes 1961 
#                        4 Liauw and Kwan 1984
#                        5 Mainstone 1974 (Adopted in FEMA 306)
#                        6 Stafford  Smith 1961 : there are 2 different options (one is equal to case 2 Paulay and Priestley 1992) and the other is 0.1*dw 
#                        7 Decanini-Fantin uncracked masonry infill
#                        8 Decanini-Fantin cracked masonry infill
# 1 Bertoldi 1993
if {$StrutWidthflag==1} { 
if {$lambdaH < 3.14} {set K1 1.300; set K2 -0.178}
if {$lambdaH > 3.14 && $lambdaH < 7.85} {set K1 0.707; set K2  0.010}
if {$lambdaH > 7.85} {set K1 0.470; set K2  0.040}
set bw		[expr $dw*($K1/$lambdaH+$K2)*$rp]; #mm , equivalent width of the strut, where rp was introduced to considered the case of opening (reduction of strut width) 
set Aw		[expr $bw*$tw*$mm*$mm]; #m2 , Area of equivalent strut
set ratio   [expr $bw/$dw];
puts [format "Strut parameter Bertoldi 1993 :k1=%.3f k2=%.3f " $K1 $K2 ]
puts [format "Strut parameter Bertoldi 1993 :bw mm =%.3f Aw mm2=%.3f bw/dw=%.3f " $bw $Aw $ratio]
}
# ----------------------------------------------------------------------------------
 # 2 Paulay and Priestley 1992
if {$StrutWidthflag==2} {
set bw      [expr $dw*0.25*$rp]; #mm , equivalent width of the strut
set Aw		[expr $bw*$tw*$mm*$mm]; #m2 , Area of equivalent strut
set ratio   [expr $bw/$dw];

puts [format "Strut parameter Paulay and Priestley 1992: bw=%.3f Aw=%.3fbw/dw=%.3f " $bw $Aw $ratio]
}
# ----------------------------------------------------------------------------------
# 3 Holmes 1961
if {$StrutWidthflag==3} {
set bw      [expr $dw*0.33*$rp]; #mm , equivalent width of the strut 
set Aw	    [expr $bw*$tw*$mm*$mm]; #m2 , Area of equivalent strut
set ratio   [expr $bw/$dw];

puts [format "Strut parameter Holmes 1961: bw=%.3f Aw=%.3f bw/dw=%.3f " $bw $Aw $ratio]
}
# ----------------------------------------------------------------------------------
if {$StrutWidthflag==4} {
set bw      [expr 0.95*sin(2*$theta)/2/sqrt($lambdaH)*$dw*$rp];
set Aw      [expr $bw*$tw*$mm*$mm]; #m2 , Area of equivalent strut
set ratio   [expr $bw/$dw];

puts [format "Strut parameter Liauw and Kwan 1984: bw=%.3f Aw=%.3f bw/dw=%.3f " $bw $Aw $ratio]
}
# ----------------------------------------------------------------------------------
# 5 Mainstone 1974 (Adopted in FEMA 306)
if {$StrutWidthflag==5} {
set bw      [expr 0.175*pow($lambdaH,-0.4)*$dw*$rp]; #mm , equivalent width of the strut
set Aw      [expr $bw*$tw*$mm*$mm];              #m2 , Area of equivalent strut
set ratio   [expr $bw/$dw];

puts [format "Strut parameter Mainstone 1974 Adopted in FEMA 306: bw=%.3f Aw=%.3f bw/dw=%.3f " $bw $Aw $ratio]
}
# ----------------------------------------------------------------------------------
# 6 Stafford  Smith 1961
if {$StrutWidthflag==6} {
set bw      [expr 0.1*$dw*$rp]; #mm , equivalent width of the strut
set Aw		[expr $bw*$tw*$mm*$mm]; #m2 , Area of equivalent strut
set ratio   [expr $bw/$dw];

puts [format "Strut parameter Stafford  Smith 1961: bw=%.3f Aw=%.3f bw/dw=%.3f " $bw $Aw $ratio]
}
# ----------------------------------------------------------------------------------
# 7 Decanini-Fantin Uncracked masonry infill
if {$StrutWidthflag==7} {
if {$lambdaH <= 7.85} {
set bw      [expr (0.085+0.748/$lambdaH)*$dw*$rp]; #mm , equivalent width of the strut
} else {
set bw      [expr (0.130+0.393/$lambdaH)*$dw*$rp]; #mm , equivalent width of the strut
}
set Aw		[expr $bw*$tw*$mm*$mm]; #m2 , Area of equivalent strut
set ratio   [expr $bw/$dw];

puts [format "Decanini-Fantin Uncracked masonry infill: bw=%.3f Aw=%.3f bw/dw=%.3f " $bw $Aw $ratio]
}
# ----------------------------------------------------------------------------------
# 8 Decanini-Fantin Cracked masonry infill
if {$StrutWidthflag==8} {
if {$lambdaH <=7.85} {
set bw      [expr (0.010+0.707/$lambdaH)*$dw*$rp]; #mm , equivalent width of the strut
} else {
set bw      [expr (0.040+0.470/$lambdaH)*$dw*$rp]; #mm , equivalent width of the strut
}
set Aw		[expr $bw*$tw*$mm*$mm]; #m2 , Area of equivalent strut
set ratio   [expr $bw/$dw];

puts [format "Decanini-Fantin Cracked masonry infill: bw=%.3f Aw=%.3f bw/dw=%.3f " $bw $Aw $ratio]
}
# ----------------------------------------------------------------------------------
# 9 Cavalieri et al.
if {$StrutWidthflag==9} {
set a1		 [expr pow(sin($theta),4)+pow(cos($theta),4)]
set a2		 [expr pow(sin($theta)*cos($theta),2)]
set a3		 [expr (pow($Ewh,-1) + pow($Ewv,-1)- pow($Gw,-1))*$a2]
set a4       [expr $v/$Ewh*$a1]
set a5       [expr $a4-$a3]
set lambda   [expr $Ewtheta/$Ec*$tw*$H/$Ac*((pow($H,2)/pow($B,2)) + (0.25*$Ac*$B/$Ab/$H))]
set vd       [expr $a5*$Ewtheta]
set c        [expr 0.249-0.0116*$vd+0.567*pow($vd,2)];
set beta	 [expr 0.146+0.0073*$vd+0.126*pow($vd,2)];
set zeta	 [expr 0.250*$B/$H -0.250+1.00]
#set epsilonv [expr $Fv/2/$Ac/$Ec]; needs to define as input the axial force acting on the column otherwise is 0
set k 1.00 ; # must be modified if is the aim is to take into account the effect of the loads k=1+(18*$lambdaH+200)*$epsilonv
set bw	    [expr $k*$c*$dw/$zeta/pow($lambda,$beta)*$rp]
set Aw		[expr $bw*$tw*$mm*$mm]; #m2 , Area of equivalent strut
set ratio   [expr $bw/$dw];
puts [format "Cavalieri et al.: lambda=%.3f" $lambda]
puts [format "Cavalieri et al.: bw=%.3f Aw=%.3f bw/dw=%.3f " $bw $Aw $ratio]
}
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# DETERMINE FAILURE MECHANISM USING DIFFERENT ANALITICAL RELATIONSHIP 
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
#put flag option to Set: 
#        FAILURE MECHANISM: $Criticalstressflag==
#                        1 Decanini Fantini-Bertoldi et al. 1987
#                        2 Decanini Fantini-Bertoldi et al. 1987 with a reduction of fwu (shear sliding of bed joint ) according to Morandi et al. and 2018 and EC6
#                        3 Paulay Priestley and Calvi 
#                        4 Paulay Priestley and Calvi  with a reduction of fwu (shear sliding of bed joint ) according to Morandi et al. 2018. and 2018 and EC6
#                        5 FEMA 306
#                        6 Eurocode EC8-Part 1/EC6 but without Fvo REDUCTION
#                        7 Eurocode EC8-Part 1/EC6 with Fvo REDUCTION
# -------------------------------------------------------------------------------
# Decanini Fantini-Bertoldi et al. 1987
if {$Criticalstressflag==1} {
if {$lambdaH < 3.14} {set K1 1.300; set K2 -0.178}
if {$lambdaH > 3.14 && $lambdaH < 7.85} {set K1 0.707; set K2  0.010}
if {$lambdaH > 7.85} {set K1 0.470; set K2  0.040}
set sigw1 [expr 1.16*$fwv*tan($theta)/($K1+$K2*$lambdaH)]; # Compression in centre
set sigw2 [expr 1.12*$fwv*sin($theta)*cos($theta)/($K1*pow($lambdaH,-0.12)+$K2*pow($lambdaH,0.88))]; # Compression at corners
set sigw3 [expr ($fwu*(1.2*sin($theta)+0.45*cos($theta))+0.3*$sig_v)*$dw/$bw]; # Shear sliding
if {$fws>0} {
set sigw4 [expr (0.6*$fws+0.3*$sig_v)*$dw/$bw]; # Diagonal tension
}
#Failure mechanism 
if {$fws>0} {
set sigw [expr min($sigw1,$sigw2,$sigw3,$sigw4)]
} else { 
set sigw        [expr min($sigw1,$sigw2,$sigw3)] 
}
set sigma		[expr $sigw*$MPa];
set force	    [expr $sigma*$Aw];

puts [format "Failure mechanism Decanini Fantin"]
puts [format "Compression in centre MPa=%.3f" $sigw1]
puts [format "Compression at corners MPa= %.3f" $sigw2]
puts [format "Shear sliding= %.3f" $sigw3]
if {$fws>0} {
puts [format "Diagonal tension= %.3f" $sigw4]
}
puts [format "Maximum stress tollerable MPa = %.3f" $sigma]
puts [format "Maximum Force tollerable kN= %.3f" $force]
set Ho [expr $force*cos($theta)];
puts [format "Horizontal shear force: "]
puts [format "Ho= %.3f " $Ho]
}
# -------------------------------------------------------------------------------
# Decanini Fantini with a reduction of fwu (shear sliding of bed joint ) according to Morandi et al. 2018 and 2018 and EC6
if {$Criticalstressflag==2} {
set fwu [expr $fwu/2];
if {$lambdaH < 3.14} {set K1 1.300; set K2 -0.178}
if {$lambdaH > 3.14 && $lambdaH < 7.85} {set K1 0.707; set K2  0.010}
if {$lambdaH > 7.85} {set K1 0.470; set K2  0.040}
set sigw1 [expr 1.16*$fwv*tan($theta)/($K1+$K2*$lambdaH)]; # Compression in centre
set sigw2 [expr 1.12*$fwv*sin($theta)*cos($theta)/($K1*pow($lambdaH,-0.12)+$K2*pow($lambdaH,0.88))]; # Compression at corners
set sigw3 [expr ($fwu*(1.2*sin($theta)+0.45*cos($theta))+0.3*$sig_v)*$dw/$bw]; # Shear sliding
if {$fws>0} {
set sigw4 [expr (0.6*$fws+0.3*$sig_v)*$dw/$bw]; # Diagonal tension
}
#Failure mechanism 
if {$fws>0} {
set sigw [expr min($sigw1,$sigw2,$sigw3,$sigw4)]
} else { 
set sigw        [expr min($sigw1,$sigw2,$sigw3)] 
}
set sigma		[expr $sigw*$MPa];
set force	    [expr $sigma*$Aw];
puts [format "Failure mechanism Decanini Fantin with reduction of fwu"]
puts [format "Compression in centre MPa=%.3f" $sigw1]
puts [format "Compression at corners MPa= %.3f" $sigw2]
puts [format "Shear sliding= %.3f" $sigw3]
if {$fws>0} {
puts [format "Diagonal tension= %.3f" $sigw4]
}
puts [format "Maximum stress tollerable MPa = %.3f" $sigma]
puts [format "Maximum Force tollerable kN= %.3f" $force]
set Ho [expr $force*cos($theta)];
puts [format "Horizontal shear force: "]
puts [format "Ho= %.3f " $Ho]
}
# ---------------------------------------------------------------------------------------------
#Paulay Priestley and Calvi 
if {$Criticalstressflag==3} {
set Vs [expr $fwu/(1-0.3*$H/$B)*$Bw*$tw/1000]; # Shear sliding failure of the infill [kN]
set Vc [expr 2*$z*$tw*$flat/3/1000]; # Compression failure of diagonal strut[kN]
set Vt [expr $pi/2*$tw*$dw*$fwv/1000]; # diagonal tension cracking [kN]
set force [expr min($Vs,$Vc,$Vt)];

puts [format "Failure mechanism Paulay Priestley and Calvi:"]
puts [format "Shear sliding=%.3f" $Vs] 
puts [format "Diagonal compression= %.3f" $Vc]
puts [format "Diagonal tension= %.3f" $Vt] 
puts [format "Force= %.3f " $force]
set Ho [expr $force*cos($theta)];
puts [format "Horizontal shear force: "]
puts [format "Ho= %.3f " $Ho]
}
# ---------------------------------------------------------------------------------------------
#Paulay Priestley and Calvi with a reduction of fwu (shear sliding of bed joint ) according to Morandi et al. 2018 and 2018 and EC6
if {$Criticalstressflag==4} {   
set fwu [expr ($fwu/2)];
set Vs [expr $fwu/(1-0.3*$H/$B)*$Bw*$tw/1000]; # Shear sliding failure of the infill [kN]
set Vc [expr 2*$z*$tw*$flat/3/1000]; # Compression failure of diagonal strut[kN]
set Vt [expr $pi/2*$tw*$dw*$fwv/1000]; # diagonal tension cracking [kN]
set force [expr min($Vs,$Vc,$Vt)];

puts [format "Failure mechanism Paulay Priestley and Calvi:"]
puts [format "Shear sliding=%.3f" $Vs] 
puts [format "Diagonal compression= %.3f" $Vc]
puts [format "Diagonal tension= %.3f" $Vt] 
puts [format "Force= %.3f " $force]
set Ho [expr $force*cos($theta)];
puts [format "Horizontal shear force: "]
puts [format "Ho= %.3f " $Ho]
}
# ---------------------------------------------------------------------------------------------
#FEMA 306
if {$Criticalstressflag==5} {   
set bw   [expr 0.175*pow($lambdaH,-0.4)*$dw]; #mm , equivalent width of the strut
set Vs   [expr ($fwu+(0.4*$sig_v))*$Bw*$tw/1000]; # Shear sliding failure of the infill [kN]
set Vc   [expr $bw*$tw*$flat*cos($theta)/1000]; # Compression failure [kN]
if {$sigm_cr>0} {
set Vcr  [expr 2*pow(2,0.5)*$tw*$Bw*$sigm_cr/1000*$Hw*$Bw/(pow($Bw,2)+pow($Hw,2))];# Diagonal cracking failure [kN]
}
set Vmi  [expr 2*pow(0.0069*$fwv,0.5)*$tw*$Bw/1000];# [kN]
set Vmf  [expr 0.3*$Vmi];# [kN]
if {$sigm_cr>0} {
set force [expr min($Vs,$Vc,$Vcr)];
} else {
set force [expr min($Vs,$Vc)];
}
if {$force>=$Vmi} {
set force $Vmi
}
if {$force<=$Vmf} {
set force $Vmf
}

puts [format "Failure mechanism FEMA 306:"]
puts [format "Shear sliding=%.3f" $Vs] 
puts [format "Compression failure= %.3f" $Vc]
if {$sigm_cr>0} {
puts [format "Diagonal cracking failure= %.3f" $Vcr] 
}
puts [format "Lim. Max strength= %.3f" $Vmi]
puts [format "Lim. Min strength= %.3f" $Vmf] 
puts [format "Force= %.3f " $force]
set Ho [expr $force*cos($theta)];
puts [format "Horizontal shear force: "]
puts [format "Ho= %.3f " $Ho]
}
# ---------------------------------------------------------------------------------------------
#Eurocode EC8-Part 1/EC6 but without Fvo REDUCTION
if {$Criticalstressflag==6} {
set force [expr ($fwu+0.4*$sig_v)*$Bw*$tw/1000]; # Shear sliding failure of the infill [kN]

puts [format "Failure mechanism Eurocode EC8-Part 1 : "]
puts [format "Force= %.3f " $force]
set Ho [expr $force*cos($theta)];
puts [format "Horizontal shear force: "]
puts [format "Ho= %.3f " $Ho]
}

#Eurocode EC8-Part 1/EC6 with Fvo reduction
if {$Criticalstressflag==7} {
set force [expr (0.5*$fwu+0.4*$sig_v)*$Bw*$tw/1000]; # Shear sliding failure of the infill [kN]

puts [format "Failure mechanism Eurocode EC6 with Fvo reduction: "]
puts [format "Force= %.3f " $force]
set Ho [expr $force*cos($theta)];
puts [format "Horizontal shear force: "]
puts [format "Ho= %.3f " $Ho]
}
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# According to Guideline Reluis : Horizontal shear force 
set Ho [expr min($fwu*$Bw*$tw/0.6,0.8*$fwv*pow(cos($theta),2)*pow(($Ec/$Ewv*$Ic*$Hw*pow($tw,3)),0.25))/1000]
puts [format "Horizontal shear force according to Lineguide Reluis: "]
puts [format "Ho= %.3f " $Ho]
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------



# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# DETERMINE BACKBONE HYSTERESIS USING DIFFERENT ANALITICAL RELATIONSHIP 
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
#put flag option to Set: 
#        BACKBONE HYSTERESIS: $Backboneflag==
#                        1 Bertoldi
#                        2 Panagiotakos and Fardis
#                        3 De Risi et al. 
#                        4 Bertoldi + Sassun et. al 
# ------------------------------------------------------------------------------
#Bertoldi
if {$Backboneflag==1} {   
set alphapost 0.02 ; 
#Stiffness
set Ksec   [expr $Ewtheta*$bw*$tw/$dw*pow(cos($theta),2)];
set Kel    [expr 4*$Ksec];
set Kpost  [expr $alphapost*$Kel];
#Force
set force1   [expr $force*0.8];
set force2   [expr $force*1.1];
set force3   [expr $force2*0.35];
#Displacement
set u_cr     [expr $force1/$Kel];
set u_max    [expr $force2/$Ksec];
set u_res    [expr $u_max + (-$force3+$force2)/$Kpost];
puts [format "Backbone curve Bertoldi: 
Stifness 
Kel=%.3f  Ksec= %.3f  Kpost= %.3f" $Kel $Ksec $Kpost]
puts [format "Force kN: 
Force1= %.3f Force2= %.3f Force3= %.3f" $force1 $force2 $force3]
puts [format "Displacement: 
Displacement cr= %.3f Displacement max= %.3f Displacement residual= %.3f" $u_cr $u_max $u_res]
#set u_zero 0.0
#set f_zero 0.0
#set i 1
#set results [list $u_zero $u_cr $u_max $u_res $f_zero $force1 $force2 $force3]
#set f [open "r$i.txt w"]
#puts $f [format $results]
#close $f
set sigma1 [expr $force1/$Aw]
set sigma2 [expr $force2/$Aw]
set sigma4 [expr $force3/$Aw]
set dw    [expr $dw*$mm]
set epsilon1 [expr $u_cr/$dw]
set epsilon2 [expr $u_max/$dw]
set epsilon3 [expr $u_res/$dw]
set epsilon4 [expr $u_res/$dw]
#Positive envelope, tension behaviour of masonry infill----- Tension/Compression==1/10

set sigma1_p [expr $force1/$Aw/10]
set sigma2_p [expr $force2/$Aw/10]
set sigma4_p [expr $force3/$Aw/10]

set dw         [expr $dw*$mm]
set epsilon1_p [expr $u_cr/$dw/10]
set epsilon2_p [expr $u_max/$dw/10]
set epsilon4_p [expr $u_res/$dw/10]


puts [format "Sigma1 kN/m2= %.3f" $sigma1]
puts [format "Sigma2  kN/m2= %.3f" $sigma2]
puts [format "Sigma4  kN/m2= %.3f" $sigma4]
}
# ------------------------------------------------------------------------------
#Panagiotakos and Fardis 
if {$Backboneflag==2} {
set alphapost 0.01 ;
#Stiffness
set Kel    [expr $Gw*$Bw*$tw/$Hw];# controllare se effettivamente la relazione per Kel è questa --.... Kel=Ew*te*bw/dw
set Kpost     [expr $Ewtheta*$bw*$tw/$dw];
set Ksoft  [expr $alphapost*$Kel];
#Force
set force1   $force;
set force2   [expr $force1*(1.3)];
set force3   [expr $force2*0.1];
#Displacement
set u_cr     [expr $force1/$Kel];
set u_max    [expr $u_cr + ($force2-$force1)/$Kpost];
set u_res    [expr $u_max + (-$force3+$force2)/$Ksoft];
puts [format "Backbone curve Panagiotakos and Fardis : 
Stifness 
Kel=%.3f  Kpost= %.3f  Ksoft= %.3f" $Kel $Kpost $Ksoft]
puts [format "Force kN: 
Force1= %.3f Force2= %.3f Force3= %.3f" $force1 $force2 $force3]
puts [format "Displacement: 
Displacement cr= %.3f Displacement max= %.3f Displacement residual= %.3f" $u_cr $u_max $u_res]
set u_zero 0.0
set f_zero 0.0
#set results [list $u_zero $u_cr $u_max $u_res $f_zero $force1 $force2 $force3]
#set f [open result.txt w]
#puts $f [format $results]
#close $f
set sigma1 [expr $force1/$Aw]
set sigma2 [expr $force2/$Aw]
set sigma4 [expr $force3/$Aw]
set dw    [expr $dw*$mm]
set epsilon1 [expr $u_cr/$dw]
set epsilon2 [expr $u_max/$dw]
set epsilon3 [expr $u_res/$dw]
set epsilon4 [expr $u_res/$dw]
#Positive envelope, tension behaviour of masonry infill----- Tension/Compression==1/10

set sigma1_p [expr $force1/$Aw/10]
set sigma2_p [expr $force2/$Aw/10]
set sigma4_p [expr $force3/$Aw/10]

set dw         [expr $dw*$mm]
set epsilon1_p [expr $u_cr/$dw/10]
set epsilon2_p [expr $u_max/$dw/10]
set epsilon4_p [expr $u_res/$dw/10]


puts [format "Sigma1 kN/m2= %.3f" $sigma1]
puts [format "Sigma2  kN/m2= %.3f" $sigma2]
puts [format "Sigma4  kN/m2= %.3f" $sigma4]
}
# ------------------------------------------------------------------------------
#Panagiotakos and Fardis Modified by De Risi et al.
if {$Backboneflag==3} {
set alphapost 0.1 ;
#Stiffness
set Ksec   [expr 0.8*$Ewtheta*$bw*$tw/$dw*pow(cos($theta),2)];
set Kel    [expr $Ksec*2.8]
set Kpost  [expr $alphapost*$Ksec];
#Force
set force1   [expr $force*0.7];
set force2   $force
set force3   [expr $force2*0.00];
#Displacement
set u_cr     [expr $force1/$Kel];
set u_max    [expr $force2/$Ksec];
set u_res    [expr $u_max + (-$force3+$force2)/$Kpost];
puts [format "Backbone curve De Risi et al.: 
Stifness 
Kel=%.3f  Ksec= %.3f  Kpost= %.3f" $Kel $Ksec $Kpost]
puts [format "Force kN: 
Force1= %.3f Force2= %.3f Force3= %.3f" $force1 $force2 $force3]
puts [format "Displacement: 
Displacement cr= %.3f Displacement max= %.3f Displacement residual= %.3f" $u_cr $u_max $u_res]
set u_zero 0.0
set f_zero 0.0
#set results [list $u_zero $u_cr $u_max $u_res $f_zero $force1 $force2 $force3]
#set f [open result.txt w]
#puts $f [format $results]
#close $f
set sigma1 [expr $force1/$Aw]
set sigma2 [expr $force2/$Aw]
set sigma4 [expr $force3/$Aw]
set dw    [expr $dw*$mm]
set epsilon1 [expr $u_cr/$dw]
set epsilon2 [expr $u_max/$dw]
set epsilon3 [expr $u_res/$dw]
set epsilon4 [expr $u_res/$dw]
#Positive envelope, tension behaviour of masonry infill----- Tension/Compression==1/10

set sigma1_p [expr $force1/$Aw/10]
set sigma2_p [expr $force2/$Aw/10]
set sigma4_p [expr $force3/$Aw/10]

set dw         [expr $dw*$mm]
set epsilon1_p [expr $u_cr/$dw/10]
set epsilon2_p [expr $u_max/$dw/10]
set epsilon4_p [expr $u_res/$dw/10]


puts [format "Sigma1 kN/m2= %.3f" $sigma1]
puts [format "Sigma2  kN/m2= %.3f" $sigma2]
puts [format "Sigma4  kN/m2= %.3f" $sigma4]
}
# ------------------------------------------------------------------------------
#Bertoldi + Limite state and performance criteria based on Sassun et al. 2015
# ------------------------------------------------------------------------------
if {$Backboneflag==4} {
set alphapost 0.02 ; 
#Stiffness
set Ksec   [expr $Ewtheta*$bw*$tw/$dw*pow(cos($theta),2)];
set Kel    [expr 4*$Ksec];
set Kpost  [expr $alphapost*$Kel];
#Force
set force1   [expr $force*1];
set force2   [expr $force*1.25];
set force3   [expr $force2*0.25];
set thetaDS1	0.0018
set thetaDS2	0.0046
set thetaDS3	0.0105
set thetaDS4 	0.0188

set u_cr	[expr $thetaDS1*$H/1000/2];
set u_max 	[expr $thetaDS2*$H/1000/2];
set u_res	[expr $thetaDS4*$H/1000/2];
puts [format "Backbone curve Bertoldi-Sassun et al. 2015:
Stifness 
Kel=%.3f  Ksec= %.3f  Kpost= %.3f" $Kel $Ksec $Kpost]
puts [format "Force kN: 
Force1= %.3f Force2= %.3f Force3= %.3f" $force1 $force2 $force3]
puts [format "Displacement: 
Displacement cr= %.3f Displacement max= %.3f Displacement residual= %.3f" $u_cr $u_max $u_res]
set u_zero 0.0
set f_zero 0.0
#set results [list $u_zero $u_cr $u_max $u_res $f_zero $force1 $force2 $force3]
#set f [open result.txt w]
#puts $f [format $results]
#close $f
set sigma1 [expr $force1/$Aw]
set sigma2 [expr $force2/$Aw]

set sigma4 [expr $force3/$Aw]
set dw    [expr $dw*$mm]

set epsilon1 0.0006
set epsilon2 0.0013
set epsilon3 0.0045
set epsilon4 0.0045
#Positive envelope, tension behaviour of masonry infill----- Tension/Compression==1/10

set sigma1_p [expr $force1/$Aw/10]
set sigma2_p [expr $force2/$Aw/10]
set sigma4_p [expr $force3/$Aw/10]

set dw         [expr $dw*$mm]
set epsilon1_p [expr $u_cr/$dw/10]
set epsilon2_p [expr $u_max/$dw/10]
set epsilon4_p [expr $u_res/$dw/10]

puts [format "Sigma1 kN/m2= %.3f" $sigma1]
puts [format "Sigma2  kN/m2= %.3f" $sigma2]
puts [format "Sigma4  kN/m2= %.3f" $sigma4]
}


# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# DETERMINE THE  HYSTERETIC BEHAVIOUR OF MASONRY INFILL USING 2 DIFFERENT MODEL
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
#put flag option to Set: 
#        HYSTERETIC BEHAVIOUR OF MASONRY INFILL: $Hysteresticflag==
#                        1 Coefficients for Unloading/Reloading Stiffness degradation and Coefficients for Strength degradation ==0
#                        2 Coefficients for Unloading/Reloading Stiffness degradation and Coefficients for Strength degradation based on N. Mohammad Noh et al. / Engineering Structures 150 (2017) 599–621
# ------------------------------------------------------------------------------
#Coefficients for Unloading/Reloading Stiffness degradation and Coefficients for Strength degradation ==0

if {$Hysteresticflag==1} {
set pF [list $sigma1_p $sigma2_p $sigma4_p $sigma4_p];
set nF [list -$sigma1 -$sigma2 -$sigma4 -$sigma4];
set pD [list $epsilon1_p $epsilon2_p $epsilon4_p $epsilon4_p];
set nD [list -$epsilon1 -$epsilon2 -$epsilon3 -$epsilon4];
set rDisp [list 0.8 0.8]; # Pos_env. Neg_env.##### Ratio of maximum deformation at which reloading begins
set rForce [list 0.1 0.1]; # Pos_env. Neg_env.##### Ratio of envelope force (corresponding to maximum deformation) at which reloading begins
set uForce [list 0.0 0.0]; # Pos_env. Neg_env.##### Ratio of monotonic strength developed upon unloading
set gammaK [list 0.0 0.0 0.0 0.0 0.0]; # gammaK1 gammaK2 gammaK3 gammaK4 gammaKLimit
set gammaD [list 0.0 0.0 0.0 0.0 0.0]; # gammaD1 gammaD2 gammaD3 gammaD4 gammaDLimit
set gammaF [list 0.0 0.0 0.0 0.0 0.0]; # gammaF1 gammaF2 gammaF3 gammaF4 gammaFLimit
set gammaE 0.0
set damage "energy"; # damage type (option: "energy", "cycle")
}
# ------------------------------------------------------------------------------
#Coefficients for Unloading/Reloading Stiffness degradation and Coefficients for Strength degradation based on N. Mohammad Noh et al. / Engineering Structures 150 (2017) 599–621
if {$Hysteresticflag==2} {
set pF [list $sigma1_p $sigma2_p $sigma4_p $sigma4_p];
set nF [list -$sigma1 -$sigma2 -$sigma4 -$sigma4];
set pD [list -$epsilon1_p -$epsilon2_p -$epsilon4_p -$epsilon4_p];
set nD [list -$epsilon1 -$epsilon2 -$epsilon3 -$epsilon4];
set rDisp [list 0.8 0.8]; # Pos_env. Neg_env. ##### Ratio of maximum deformation at which reloading begins
set rForce [list 0.1 0.1]; # Pos_env. Neg_env. ##### Ratio of envelope force (corresponding to maximum deformation) at which reloading begins
set uForce [list 0.0 0.0]; # Pos_env. Neg_env. ##### Ratio of monotonic strength developed upon unloading
set gammaK [list 0.8 0.7 0.7 0.7 0.8]; # gammaK1 gammaK2 gammaK3 gammaK4 gammaKLimit ##### Coefficients for Unloading Stiffness degradation
set gammaD [list 0.0 0.0 0.0 0.0 0.1]; # gammaD1 gammaD2 gammaD3 gammaD4 gammaDLimit ##### Coefficients for Reloading Stiffness degradation
set gammaF [list 1 0.0 1 1 0.1]; # gammaF1 gammaF2 gammaF3 gammaF4 gammaFLimit ##### Coefficients for Strength degradation
set gammaE 10
set damage "cycle"; # damage type (option: "energy", "cycle")
}

# add the material to domain through the use of a procedure
procUniaxialPinching 8${eleTag} $pF $nF $pD $nD $rDisp $rForce $uForce $gammaK $gammaD $gammaF $gammaE $damage
proc procUniaxialPinching {materialTag pEnvelopeStress nEnvelopeStress pEnvelopeStrain nEnvelopeStrain rDisp rForce uForce gammaK gammaD gammaF gammaE damage} {
uniaxialMaterial Pinching4 $materialTag [lindex $pEnvelopeStress 0] [lindex $pEnvelopeStrain 0] [lindex $pEnvelopeStress 1] [lindex $pEnvelopeStrain 1] [lindex $pEnvelopeStress 2] [lindex $pEnvelopeStrain 2] [lindex $pEnvelopeStress 3] [lindex $pEnvelopeStrain 3] [lindex $nEnvelopeStress 0] [lindex $nEnvelopeStrain 0] [lindex $nEnvelopeStress 1] [lindex $nEnvelopeStrain 1] [lindex $nEnvelopeStress 2] [lindex $nEnvelopeStrain 2] [lindex $nEnvelopeStress 3] [lindex $nEnvelopeStrain 3] [lindex $rDisp 0] [lindex $rForce 0] [lindex $uForce 0] [lindex $rDisp 1] [lindex $rForce 1] [lindex $uForce 1] [lindex $gammaK 0] [lindex $gammaK 1] [lindex $gammaK 2] [lindex $gammaK 3] [lindex $gammaK 4] [lindex $gammaD 0] [lindex $gammaD 1] [lindex $gammaD 2] [lindex $gammaD 3] [lindex $gammaD 4] [lindex $gammaF 0] [lindex $gammaF 1] [lindex $gammaF 2] [lindex $gammaF 3] [lindex $gammaF 4] $gammaE $damage
}

# ---------------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------------
# -------------------------ELEMENT DEFINITION ACCORDING TO OPENSEES----------------------
# ---------------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------------

if {$typ=="single"} {
set nI	[lindex $nds 0];
set nJ	[lindex $nds 1];
set nK	[lindex $nds 2];
set nL	[lindex $nds 3];
#element truss $eleTag $iNode $jNode $A $matTag <-doRayleigh $rFlag>
element truss ${eleTag}1 $nI $nK $Aw 8${eleTag} -doRayleigh 1
element truss ${eleTag}2 $nL $nJ $Aw 8${eleTag} -doRayleigh 1
}
puts $pfile [format "Element ${eleTag}1 Force:%.5f %.5f %.5f %.5f Displacement:%.5f %.5f %.5f %.5f " $force1 $force2 $force3 $force3 $u_cr $u_max $u_res $u_res]
}